#pragma once
#include <Geometry/cubic_bspline_curve_cuda.h>
#include <device_launch_parameters.h>
namespace shermit {
namespace fitting {
namespace curve {

	__global__ void calcBaseFuncMatrix32(
		CubicBSplineCurve<float, DeviceType::GPU> curve,
		int nd, int nc,
		const float* t,
		float* matrix
	)
	{
		int idx = blockIdx.x * blockDim.x + threadIdx.x;
		if (idx >= nd * nc) return;
		int i = idx / nd;      // column index
		int j = idx - i * nd;  // row index
		matrix[idx] = curve.baseFuncVal(i, t[j]);
	}
	__global__ void calcBaseFuncMatrix64(
		CubicBSplineCurve<double, DeviceType::GPU> curve,
		int nd, int nc,
		const double* t,
		double* matrix
	)
	{
		int idx = blockIdx.x * blockDim.x + threadIdx.x;
		if (idx >= nd * nc) return;
		int i = idx / nd;      // column index
		int j = idx - i * nd;  // row index
		matrix[idx] = curve.baseFuncVal(i, t[j]);
	}

	__global__ void calcBaseFuncMatrixClosed64(
		CubicBSplineCurve<double, DeviceType::GPU> curve,
		int nd, int n,
		const double* t,
		double* matrix
	)
	{
		int idx = blockIdx.x * blockDim.x + threadIdx.x;
		if (idx >= nd * n) return;
		int i = idx / nd;      // column index
		int j = idx - i * nd;  // row index
		matrix[idx] = curve.baseFuncVal(i, t[j]);
		if (i < 3) {
			matrix[idx] += curve.baseFuncVal(n + i, t[j]);
		}
	}

	__global__ void calcPtSOA32(
		int size,
		const float *base_func_mtx, 
		const nv::vec2<float>* cpts,
		const nv::vec2<float>* pts, float* pt_soa) 
	{
		int i = blockIdx.x * blockDim.x + threadIdx.x;
		if (i >= size) return;
		pt_soa[i] = pts[i].x;
		pt_soa[i + size] = pts[i].y;
	}

	
	__global__ void calcPtSOA64(
		int size,
		const nv::vec2<double>* pts, double* pt_soa)
	{
		int i = blockIdx.x * blockDim.x + threadIdx.x;
		if (i >= size) return;
		pt_soa[i] = pts[i].x;
		pt_soa[i + size] = pts[i].y;
	}

	__global__ void setPts32(int size, const float* ptx, const float* pty, nv::vec2<float>* pts) {
		int i = blockIdx.x * blockDim.x + threadIdx.x;
		if (i >= size) return;
		pts[i].x = ptx[i];
		pts[i].y = pty[i];
	}
	__global__ void setPts64(int size_out, int size_in, const double* ptx, const double* pty, nv::vec2<double>* pts) {
		int io = blockIdx.x * blockDim.x + threadIdx.x;
		if (io >= size_out) return;
		int ii = io % size_in;
		double x = ptx[ii];
		double y = pty[ii];
		if (isfinite(x) && isfinite(y)) {
			pts[io].x = x;
			pts[io].y = y;
		}
	}

}
}
}